Doppler synchronization of pulsating phases by time delay 
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Synchronization by exchange of pulses is a widespread phenomenon, observed in flashing fireflies, 
applauding audiences and the neuronal network of the brain. Hitherto the focus has been on 
integrate-and-fire oscillators. Here we consider entirely analytic time evolution. Oscillators exchange 
narrow but finite pulses. For any non-zero time lag between the oscillators complete synchronization 
occurs for any number of oscillators arranged in interaction networks whose adjacency matrix fulfils 
some simple conditions. The time to synchronization decreases with increasing time lag. 
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INTRODUCTION 

The emergence of coherent structures in time and space 
though synchronization occurs across the entire breadth 
of science: vibrating atoms, firing neurons, flashing fire- 
flies, clapping audiences etc. and has therefore been stud- 
ied intensively from a mathematical viewpoint [1-3]. 

Synchronization is often analyzed in models which ex- 
plicitly favor phase synchronization, e.g. in the seminal 
Kuramoto model [2, 4] and in diffusively coupled mod- 
els (see e.g. [5, 6]). In these schemes the net- interaction 
between oscillators indeed vanishes in the synchronized 
state. 

However, in many cases, such as fireflies [7-9], car- 
diac cells, neuronal system and applauding audiences 
[3, 10] the interaction between oscillators consists in the 
exchange of brief pulses, which persist even when the 
system fully synchronize. Since Mirollo and Strogatz's 
influential 1990 paper [11] such systems are often de- 
scribed by a set of non-analytically evolving integrate- 
and-fire oscillators. Each oscillator is described by a 
load variable, which is taken to have a concave depen- 
dence on a monotonously increasing phase. When the 
load reaches a certain threshold, relaxation occurs in- 
stantaneously and a pulse is send to all connected os- 
cillators. Receiving oscillators jump discontinuously for- 
ward by a given amount. For such systems Mirollo and 
Strogatz showed [11] that full synchronization always oc- 
curs. Later Ernst, Pawelzik and Geisel [12] demonstrated 
that for excitatory-only couplings, synchronization de- 
pends on phase lag, whereas the presence of inhibitory 
couplings leads to full in-phase synchronization. 

The treatment of pulse oscillators in terms of non- 
analytic integrate-and-fire oscillators is more a tradition 
than a necessity. In the present paper we assume that 
each oscillator is represented by a phase 0i (t) whose time 



t derivative is always equal to a constant rate plus a 
sum of smooth but narrow pulses emitted by surround- 
ing oscillators coupled with strength J. Synchronization 
(asymptotically vanishing phase difference) always occurs 
for this system if pulses arrive with a non-zero time lag St 
for a very wide class of adjacencies, including the mean- 
field setting often considered in the literature. 

Given the previous results for pulse oscillators [11, 12] 
and to illuminate the more detailed discussion below it 
is natural to begin our analysis of two pulse exchanging 
phases by considering Dirac's delta pulses for the inter- 
action. 

6\{t) = Lu + jJ2^{t-St)-n)) 

6 2 {t) = Lj + jJ2S(9i{t-5t)-n)) (1) 

n£Z 

Integrating the time derivatives tells us that 0\ "jumps" 
each time O2 passes through an integer value, 9\(t) i-> 
6>i(i) + J, and vise versa for 2 - Let 0i(Q) > #2(0), it is 
straight forward to see, Fig. 1, that the two phases are 
unable to synchronize though in the case of a finite time 
delay they may leapfrog each other, as the jump of one 
oscillator can make the other skip its. 

Obviously Dirac delta pulses are unrealistic. Pulses 
emitted by real systems will have a finite width and a 
smooth time dependence (Eq. (2) below). The introduc- 
tion of smooth pulses changes the behavior in an essen- 
tial way. As will be explained below synchronization now 
takes place whenever a time lag is present, St > 0, and in 
this case complete synchronization occurs for all smooth 
pulses. 

General model - We now consider n coupled oscilla- 
tors, each described by a single degree of freedom with 
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FIG. 1: The time evolution of two oscilators (solid and dashed 
lines) exchanging pulses according to Eq. 1. Left panel: St — 
0, right panel: St — 0.5. 

i = 1, 2, . . . , n, and each with the same eigerifrcquency u: 

ei{t) = u + J2 J iM6j(t-fo)) ( 2 ) 

3 

Oscillators are coupled through an adjancency matrix 
and a feedback function a(6) which has period £. It is 
only through a (9) that periodicity is implemented: a(9i) 
describes the effect that the state of 9i (say, the flashing 
of a firefly) has on any other oscillator. As opposed to 
other models often studied in synchronization, such as 
the Kuramoto model [2], the effect of a docs not disap- 
pear in the synchronized state. 

We chose 9\ > at all times such that 9i(t) are mono- 
tonically increasing functions in time. This is achieved by 
choosing a(9) > 0. In our numerical study below we use a 
comb of normalised Gaussians with period £ and width w, 
= E^-ooexp(-(x + <) 2 /(2^ 2 )) (2™ 2 )" 1 /* = 
£ -1 7? 3 (7ra/£,exp (-2u; 2 7r 2 /£ 2 )) , the Jacobi theta func- 
tion. 

In natural systems time delay is inevitable. We show 
that 5t > is crucial for synchronization. We use 
this term in a strong sense: For any pair i,j of oscil- 
lators Hindoo 9i(t) — 9j(t) = with /iy G Z, i.e. the 
phase difference between any two oscillators converges 
to an integer multiple of the period of a, which implies 
linii^oo 9\ — 9j = 0. By inspection it is clear that for the 
synchronized state to exist indefinitely, ■ Jij — J needs 
to be independent of i, which means that if synchroniza- 
tion takes place, the difference between any 9i{t) and the 
solution 9{t) of 

9{t) =lo + Ja(9(t-5t)) (3) 

with appropriate initial conditions vanishes asymptoti- 
cally. Provided J / 0, the eigenfrequency co can be ab- 
sorbed into a, using a{9) — > a (9) + co/J. 



FIG. 2: Comparison of <j>(t) = (0i(t) - 9 2 (t))/2 from a nu- 
merical integration of Eq. (2) (filled circles) and the linear 
approximation (dashed line) Eq. (5), with <f>(t) — 0.05 and 
St = 0.01. Parameters are n = 2, Jij = 1 — Sij, oj = 2, 
^ = 1.01 and w = 0.1, so that ip « 0.3934. The inset com- 
pares of full solution and linear approximation for (f>(0) = 0.05 
and St = 0.1. 

Simple two oscillator case - We now demonstrate that 
under very general conditions the system in Eq. Eq. (2) 
will synchronize in the long time limit. First we consider 
the simple case of two oscillators, i.e. n — 2 and = 
1 — Sij. By considering 9\/92, it is easy to show that 
9\(t) — 9 2 (t) is periodic if St =0, i.e. synchronization in 
the strong sense above does not occur without time delay, 
rather, entrainment is inevitable. However, integrating 
the equation of motion Eq. (2) numerically on the basis 
of a simple Euler method suggests differently. Better 
numerical schemes, such as the Runge Kutta [13] method, 
remove the spurious synchronization, which depends on 
the integration time step and therefore hints at the role 
of the time delay effectively implemented by the forward 
derivative used in the most naive Euler scheme. 

We now analyze in detail the effect of a time delay by 
considering Eq. (2) with St > 0. A linear stability anal- 
ysis for small St and small deviations 9{ (t) — 9(t) reveals 
that any positive St leads to a synchronized state. We 
present the calculation briefly in the following for n = 2 
and Jij = 1 — Sij . 

The equations of motion of 4>(t) = (l/2)(9 1 (t) - 9 2 (t)) 
and 9 = (l/2)(9 1 (t) + 9 2 (t)) are 

kt) = \(cr(92(t-St))-c7(9 1 (t-5t))) (4a) 

kt) = ^(o-(9 2 (t-St)) + a(9 1 (t-St)))+cj (4b) 

which we study to first order in <f> and St and find (for 
details sec [14]) 

= 0(0)5^) cxp (-2St J 2 [ m d9 a ' 2 [ 6) \ (5) 
9(t) V J*M w + MO) J 
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with t(9) = Jg {0) d9(uj + 0(e))- 1 . As the integrand is 
strictly positive, synchronization takes place in this ap- 
proximation for all St > 0. Fig. 2 shows that the lin- 
earized solution is a very good approximation of the full 
system Eq. (2). 

The characteristic time to synchronization is estimated 
in the following way. Define 



i, = 



d6- 



1 



dO- 



■ + Ja(6) 

■ + Ja{6) 



i d9_ 

9 



(6a) 



(6b) 



where tp, to leading order, is the time for a(9(t)) to go 
through one period, i.e. 9(t+ip) « 9(t)+^. S corresponds 
to the integral in the exponent of Eq. (5) for 9(t) — 0(0) + 
£. As the integrand is periodic we estimate 



(0) u> 



+ Ja(9) i> 



(7) 



and rewrite Eq. (5) <f>(t) = <f>(0)(6(0)/6(t)) exp (-t/r), 
with the characteristic synchronization time 



(8) 



can 



2J 2 StS 

inversely proportional to the time delay St . 

Network of oscillators The above picture 
be extended to arbitrary coupling matrices Jy, or a 
(weighted) network adjacency matrix. The only con- 
straint is = J independent from i, similar to 
a Markov matrix. Motivated by the observation that 
the two parameters used above, 9 and <f>, are based on 
the eigenvectors of the matrix Jy = 1 — Sij studied for 
n = 2, we consider the time evolution of (i\(j)(t)), i.e. 
of "normal modes", where (i\ is the ith left eigenvec- 
tor of J, which, we assume for simplicity, has n linearly 
independent eigenvectors. For simplicity, we normalize 
(i\j) = S^. The matrix J is not necessarily symmetric so 
generally ((i|)^ 7^ \i). Due to the Markov property, there 
is a pair of left and right eigenvectors with eigenvalue J, 
which in the following is denoted by (1| and |1) = ]P" |ej) 
respectively, where \a) denotes the canonical basis of the 
R™. 

The state of the entire system is written in vector form 
as \9{i)) = Y^i \ei)9i{t). The column vector \4>{t)) is the 
deviation \<f>(t)) = \9{t)} -J(t)\l) of \9(t)) from 9(t) = 
(l\9(t)) anticipating that 9(t) represents the asymptoti- 
cally synchronised state. Following the procedure above, 
one finds 



( f m a' 2 (9') " 

x exp Ai(J - Xi)St / „;J d9' 



T{9') 



(9) 
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FIG. 3: Plot of the synchronization time estimated from the 
window averaged phase difference 6\(t) — 82(f) (average taken 
over a time period t* so that 6(t) — 9(t — t*) — £). The filled 
symbols refer to results based on Eq. (2), the empty triangles 
to Eq. (5) and the line to Eq. (8) . Parameters as in Fig. 2. 

where T = T(9(0)) and T{9) = uj + 
J (a(6) - St9{9)a'{9)^ = 9(9) + 0(St 2 ). 

The amplitudes Ai are determined by the initial pro- 
jections (i\4>(0)) = Ai. Eq. (9) also applies to i = 1, yet 
(1 )</>(£)) = by construction so that A\ = 0. The special 
case J — (so that 9 = uj to linear order) coincides with 
the limit J — > 0, where 



lim ( T ^)) \ 3 =e \M§)/u J -6t\ i v<(e) 
J^o \ To J 



(10) 



to leading order, assuming Tq = cu for simplicity. 

Since T(9) is periodic, the long-term behaviour of 
(i\4>(t)) depends crucially on the sign of the real part 
of \i(J — \i). If it is negative, the projection has an 
approximate synchronization time 



n = 



1 

UU u+Ja(§) 



-$t(\i(J - \i))St J^d9 



2 (0) 



(11) 



corresponding to Eq. (8). Here 3?(-) denotes the real part. 
The usual mean-field setup Jjj = a(l — Sij) has one eigen- 
value J = (n — l)a and n — 1 eigenvalues Ai = —a, so 
that Aj( J — Xi) = —no? has a negative real part provided 
a 2 has a positive one, i.e. in particular for all real a. The 
mean-field theory thus always synchronises, as does the 
special case of a lattice Laplacian. 

The perturbative result Eq. (9) can be compared to the 
numerical integration of the system. We used a fourth 
order Rungc-Kutta integration scheme [13, 14] and show 
in Fig. 3 that the derived synchronization time compares 
very well (for time delays up to 5% to 10% of the syn- 
chronized period) to that of the linearised result, Eq. (5) 
and to the estimate Eq. (8). 
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FIG. 4: We see that 9\ is ahead of 62, with the maximum 
of <r(02(i)) displaced by At (as indicated) to the right rela- 
tive to that of a(8i(t)), as initialised. At a given time, 62 
still has to pass through the maximum of a(02(t)) when 9i 
already has. The phase speed #1,2 (t) is essentially o"(6*2,i(i)) 
shifted by St to the right, as indicated by the dotted lines. 
As a result, the maximum of 62 nearly aligns with the max- 
imum of cr(#2(i)), i.e. 62 is fast when a (62^)) goes through 
the maximum (dashed line), thereby narrowing it. In turn, 9i 
passes very quickly through a low value of a(di(t)), relatively 
broadening in turn the maximum of cr(6i(t)). 

Mechanism - How is synchronization achieved? Fig. 4 
shows cr(0i(t)) and Qiit) as a function of t for n = 2. 
Synchronization occurs because 6*2 experiences a greater 
increase in speed by (7(6*1 (i — St)) than Q\ does by 
cr(#2(i — St)). This asymmetry comes about because 
02 is relatively fast itself when cr(02(£)) goes through its 
maximum and 9\ is relatively slow when (7(6*1 (£)) goes 
through its maximum. As a result the maximum cr(0i(t)) 
is broadened as a function of time, and (7(6*2 (t)) is nar- 
rowed (this effect is minute and thus not visible in Fig. 4). 
Therefore, the maximum of (7(6*1 it — 6)) enters into 6*2 for 
a longer time period than a fait — S)) enters into 0\, 
leading to a speedup of 6*2 relative to 6\. In summary, 
synchronization is result of oscillator i being slow or fast 
when going through the maximum of the function crfa). 
What role has the time delay in this? The time delay 
ensures that the trailing oscillator 6*2 receives a boost at 
a time when (7(6*2 it)) goes through a maximum, while the 
leading 6>i receives its boost at a time when a fait)) goes 
through a local minimum. Without the time delay, the 
effect of the speed-up and the slowdown would indeed be 
perfectly symmetric. 

We notice that the mechanism underlying the synchro- 
nization supported by Eq. (2) is a kind of Doppler effect 
that makes the received pulse change its duration when 
the sending oscillator changes its speed. 



Equation Eq. (2) provides a remarkably simple mech- 
anism for synchronization. Because oscillators lagging 
behind by a certain amount catch up in every period of 
ct(6*) by an amount of phase difference proportional to 
the phase difference at the beginning of the period, the 
model can immediately be extended to one with differ- 
ent cigenfrcquencies uii of oscillators or some variation in 
J2j Jij with i or of a and even St . An analysis of such 
extensions follows the derivation above, see [14]. It will 
generally lead to entrainmcnt. 

Synchronisation by time delay is a viable explanation 
for natural synchronisation phenomena whenever oscila- 
tors respond to the duration of the pulse received. Fire- 
flies are known to be able to change the pulse duration 
and female fireflies are sensitive to that [8]. The exact 
way a clapping audience reaches synchrony [3, 10] can 
be analyzed sufficiently accurately to establish whether 
people change the duration of the individual clap [15] in 
the process of reaching synchrony. 
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